!! Copyright (C) 2009,2010,2011,2012  Carlo de Falco
!!
!! This file is part of:
!!   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
!!
!! LDGH is free software: you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! LDGH is distributed in the hope that it will be useful, but WITHOUT
!! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
!! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
!! License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with LDGH. If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Carlo de Falco                    <carlo.defalco@gmail.com>


module mod_petscintf
  
#include <finclude/petscdef.h>

  use petscvec
  use petscmat  
  use petscksp 

  Vec, save :: x, f, xloc
  Mat, save :: J
  KSP, save :: solver
  PetscViewer, save :: viewer

  PetscErrorCode :: ierr
  PetscInt :: ii, jj
  PetscInt, save ::  mii, nii
  
contains

!-----------------------------------------------------------------------

  subroutine mod_petscintf_constructor ()
    call PetscInitialize (PETSC_NULL_CHARACTER, ierr)
  end subroutine mod_petscintf_constructor

!-----------------------------------------------------------------------

  subroutine mod_petscintf_destructor ()
    call PetscFinalize (ierr)
  end subroutine mod_petscintf_destructor

!-----------------------------------------------------------------------

  subroutine petsc_init (dim)
    
    integer, intent(in) :: dim

    !! create matrix, solution and rhs
    mii = dim
    nii = dim
  
    !! vectors
    call VecCreate (PETSC_COMM_WORLD, &
         f, ierr)

    call VecSetSizes (f, PETSC_DECIDE, &
         nii, ierr)

    call VecSetFromOptions (f, ierr)

    call VecDuplicate (f, x, ierr)

    call VecCreateSeq (PETSC_COMM_SELF, &
         nii, xloc, ierr)


    !! matrix
    call MatCreate (PETSC_COMM_WORLD, &
         J, ierr)

    call MatSetSizes (J, PETSC_DECIDE, &
         PETSC_DECIDE, nii, nii, ierr)

    call MatSetFromOptions (J, ierr) 

  end subroutine petsc_init

!-----------------------------------------------------------------------
  
  subroutine petsc_solve (dim, ir, jc, av, &
       gn, sol, rhs)

    integer, intent(in) :: &
         dim
    integer, intent(in) :: &
         ir(:), jc(:), gn(:)
    real, intent(in)    :: &
         av(:),  rhs(:)
    real, intent(inout) :: &
         sol (:)

    integer     :: &
         nz 
    PetscScalar :: &
         tmp
    PetscInt    :: &
         itmp, jtmp

    VecScatter :: &
         vs

    PetscScalar, pointer :: &
         xloc_v(:)

    nz = size (av, 1)

    call VecSet (f, 0.0, ierr)
    CHKERRQ (ierr)

    call VecSet (x, 0.0, ierr)
    CHKERRQ (ierr)

    do ii=1,dim

       itmp = gn(ii)
       tmp  = rhs(ii)

       call VecSetValue (f, itmp, tmp, &
            ADD_VALUES, ierr)
       CHKERRQ (ierr)

       tmp = sol(ii)
       call VecSetValue (x, itmp, tmp, &
            ADD_VALUES, ierr)
       CHKERRQ (ierr)

    enddo
    
    call VecAssemblyBegin (f, &
         VEC_FINAL_ASSEMBLY, ierr)
    CHKERRQ (ierr)

    call VecAssemblyEnd (f, &
         VEC_FINAL_ASSEMBLY, ierr)  
    CHKERRQ (ierr)

    call VecAssemblyBegin (x, &
         VEC_FINAL_ASSEMBLY, ierr)
    CHKERRQ (ierr)

    call VecAssemblyEnd (x, &
         VEC_FINAL_ASSEMBLY, ierr)  
    CHKERRQ (ierr)

    do ii=1,nz

       itmp = ir(ii)
       jtmp = jc(ii)
       tmp  = av(ii)

       call MatSetValue (J, itmp, jtmp, tmp, &
            ADD_VALUES, ierr)
       CHKERRQ (ierr)

    enddo

    call MatAssemblyBegin (J, &
         MAT_FINAL_ASSEMBLY, ierr)
    CHKERRQ (ierr)

    call MatAssemblyEnd (J, &
         MAT_FINAL_ASSEMBLY, ierr)
    CHKERRQ (ierr)

    !! create and apply solver
    call KSPCreate (PETSC_COMM_WORLD, &
         solver, ierr)
    CHKERRQ (ierr)

    call KSPSetOperators (solver, J, J, &
         SAME_NONZERO_PATTERN, ierr)
    CHKERRQ (ierr)

    call KSPSetFromOptions (solver, ierr)
    CHKERRQ (ierr)

    call KSPSolve (solver, f, x, ierr)
    CHKERRQ (ierr)

    !! copy solution to all processes
    call VecScatterCreateToAll (x, vs, &
         PETSC_NULL_OBJECT, ierr)
    
    call VecScatterBegin (vs, x, xloc, &
         INSERT_VALUES, SCATTER_FORWARD, ierr)
     
    call VecScatterEnd (vs, x, xloc, &
         INSERT_VALUES, SCATTER_FORWARD, ierr)
    
    call VecGetArrayF90 (xloc, xloc_v, ierr)
    do ii=1,nii
       sol(ii) = xloc_v(ii)
    end do
    call VecRestoreArrayF90 (xloc, &
         xloc_v, ierr)

    call VecScatterDestroy (vs, ierr)    
    call VecDestroy (xloc, ierr)


end subroutine petsc_solve

!-----------------------------------------------------------------------

subroutine petsc_cleanup
  
  !! free matrix, solution and rhs
  call MatDestroy (J, ierr)
  call VecDestroy (x, ierr)
  call VecDestroy (f, ierr)

  !! free solver
  call KSPDestroy(solver, ierr)

end subroutine petsc_cleanup

!-----------------------------------------------------------------------

subroutine petsc_print ()

  !! print matrix, rhs and solution
  call PetscViewerCreate(PETSC_COMM_WORLD, &
       viewer, ierr);
  call PetscViewerSetType (viewer, &
       PETSCVIEWERASCII, ierr);
  call PetscViewerSetFormat (viewer, &
       PETSC_VIEWER_ASCII_MATLAB, ierr);
  call PetscViewerFileSetName (viewer, &
       "mod_petscintf_debug.output", ierr);
  
  call MatView (J, viewer, ierr)
  CHKERRQ (ierr)
  call VecView (x, viewer, ierr)
  CHKERRQ (ierr)
  call VecView (f, viewer, ierr)
  CHKERRQ (ierr)

  call PetscViewerDestroy (viewer, ierr)

end subroutine petsc_print

!-----------------------------------------------------------------------

end module mod_petscintf

